A Configurational Bias Monte Carlo Method for Linear 

and Cyclic Peptides 

Michael W. Deem* and Joel S. Bader 
CuraGen Corporation 
322 East Main Street 
Branford, CT 06405 

February 1, 2008 



Running Title: Biased Monte Carlo of Cyclic Peptides. 



*Author to whom correspondence should be addressed. Present address: Lyman Labo- 
ratory of Physics, Harvard University, Cambridge, MA 02138. 



1 



Abstract 



In this manuscript, we describe a new configurational bias Monte Carlo technique for 
the simulation of peptides. Wc focus on the biologically relevant cases of linear and cyclic 
peptides. Our approach leads to an efficient, Boltzmann-weighted sampling of the torsional 
degrees of freedom in these biological molecules, a feat not possible with previous Monte 
Carlo and molecular dynamics methods. 
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1 Introduction 



This paper presents a new Monte Carlo method that employs biased trial moves to achieve 
an efficient sampling of the torsional degrees of freedom for linear and cyclic peptides. 

Peptides are small molecules, built from amino acids, that are of fundamental impor- 
tance in biological systems [|l|. They play key roles in signal transduction between cells, 
regulation of cell growth and differentiation, and protein localization on cell surfaces [@]. 
Peptides are thought to regulate neurotransmission, from modulating pain and thirst to 
affecting memory and emotion p, ^. They are used as a chemical defense mechanism by 
some organisms. The conus snails, for example, produce a family of highly-constrained 
peptides that include very powerful neurotoxins p. Finally, peptides are used within the 
biotechnology industry to identify antagonists blocking various abnormal enzymatic ac- 
tions or ligand- receptor interactions ||^. Cyclic or otherwise constrained peptides are often 
preferred for this application, since such molecules suffer less of a loss of configurational 
entropy upon binding A classic example is the use of the ROD peptide to block the 
GPIIb/IIIa-fibronectin interaction, reducing blood platelet aggregation p], |[. 

The properties of peptides are amenable to examination by computer experiment. An 
early study was of the alanine dipeptide, in which the potential energy surface was de- 
duced from ab initio quantum mechanical calculations [0, |TT]]. Larger peptides have been 



examined by classical simulations. Both molecular dynamics [|12[ and Monte Carlo []T3 



approaches have proven useful. The effects of the aqueous environment have been incorpo- 
rated by simple dielectric theory [0, |I6|, [1^ or by explicit inclusion of water molecules 



IS 



It has become clear, however, that the standard molecular dynamics and Monte Carlo 
methods are not capable of sampling all conformational degrees of freedom accessible at 
body temperature to the larger peptides. This problem is particularly evident for the im- 
portant case of constrained peptides. Various solutions, such as high-temperature molec- 
ular dynamics |]T^, ^ or simplified force fields |2T|], have been suggested, but these 
approaches suffer from uncontrolled approximations. A simulation method able to sample 
the relevant conformational states of peptides, particularly constrained ones, or exposed 
loops of larger proteins would be of great value. It would aid study of these molecules in bi- 
ological systems as well as facilitate structural understanding of the peptides and antibodies 
of interest to the biotechnology industry. 
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Recently, powerful Monte Carlo methods have been developed that have a greatly en- 
hanced sampling efficiency |2^, ^ |2^, These methods have been 
applied to chain molecules at low and high density Q and even at phase coexistence 
|32| , [33| , pl| . These methods all use importance sampling, or biased moves, to efficiently 



explore the free energy landscape. 

We here apply these concepts to peptide molecules. Both linear and constrained or cyclic 
peptides are treated by this method. In Sec. 2 we describe the Monte Carlo method in detail. 
Appendices describe the rigid molecular fragments from which peptides are constructed and 
provide technical details of the method. In Sec. 3 we describe the application of this method 
to the prototypical polyglycine peptides. We discuss the results in Sec. 4. The superiority 
of this method over conventional molecular dynamics and Monte Carlo is demonstrated. 
Conclusions are presented in Sec. 5. 



2 Monte Carlo Method 

We make the simplifying assumption that the intramolecular potential energy contains only 
torsional and non-bonded terms. That is, bond lengths and angles are fixed, and rotation 
is allowed only about sigma bonds. At room- or body-temperature, these are fairly good 
assumptions. They could easily be relaxed, although sampling the increased degrees of 
freedom would entail a computational expense. Appendix A describes the rigid fragments 
that occur in peptides under these assumptions. A suitable form for the interatomic po- 



tential would be the AMBER |3|, ECCEP [||], or CHARMm force field. We pick 
the AMBER potentials. Water is treated in an implicit way, assuming the dielectric con- 
stant for Coulomb interactions is given by e/eo = 4r, with r given in Angstroms. These 
assumptions allow the method to be presented without a discussion of detailed force field 
issues. The method is generically applicable to better force fields and an explicit treatment 
of water. 

A configurational bias Monte Carlo (CBMC) technique is used to explore the conforma- 
tions of the molecules. We describe the algorithm for both linear and cyclic peptides. By 
cyclic, we mean peptides constrained because of disulfide bonds between cystine residues. 

There are two types of atoms in a peptide, those in the side chains and those in the 
backbone. Consequently, there are two types of Monte Carlo moves: type I moves change 
the positions of side chain atoms only, and type II moves change the positions of backbone 
atoms, rigidly rotating the attached side chains. The type I move is an extension of the 
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chain-molecule CBMC ^ to the structurally more complicated case of peptides. The 
type I move is applicable to side chains with a free end {i.e. all naturally occurring amino 
acid side chains except for proline). The backbone to which the side chain is attached can 
be either linear or cyclic. In the cyclic case, the type I move is also used to change the 
configuration of the free ends of the main chain. 

There are two kinds of type II moves for the backbone: type Ila moves for linear peptides 
and type lib moves for cyclic peptides. The type Ila move is essentially the same as a type 
I move. The side-chain residues that are attached to the backbone are rigidly rotated so 
as to remain properly bonded to the Ca atoms in their new positions. When the peptide 
is cyclic, we use a type lib move to change the configuration of part of the backbone loop, 
rigidly rotating any side chains or free ends of the peptide that are attached to that part of 
the backbone. The backbone of a cyclic peptide includes the atoms along the main chain 
as well as the C/3 and S atoms of the cystines participating in the disulfide bond. This 
move requires a concerted rotation of the backbone torsional angles with a rigid rotation 
of the attached side groups. This concerted rotation of the torsional angles is an extension 



of the concerted rotation scheme for alkanes |22, 28 



A type I move is initiated by identifying the side chain to be regrown. Not all of the side 
chain need be regrown, and the first group to regrow is chosen. This feature is helpful for 
the amino acids with longer side chains, such as lysine. These choices are made randomly. 
The M rigid units to be regrown are first removed and then added one at a time, starting 
from the one closest to the backbone. For each addition, the following actions are carried 
out (see Fig. 1): 

1) k values of the torsional angle (pij, 1 < j < k connecting rigid unit i to unit i — 1 are 
generated according to the internal potential, 

p^\<|>^,) cx exp[-/3«r*(0..)] • (1) 

The function u^j^^{(f)ij) is the part of the internal energy that couples unit i to the rest of the 
molecule (but excluding units i + 1 to M). The inverse temperature is given hj f3 = l/ksT. 

2) One of these is picked with probability 

pr*(0^,) = exp[-(3u^\<P^,)]/w^''\^) , (2) 

where 

k 

ext ( ■\ \ ^ fin,^^^ 



W [I 



eM-^^uT\<P^,)] . (3) 
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The function u'^^'^{(j)ij) is the part of the external energy that couples unit i to the rest of 
the molecule (but excluding units i + 1 to M). 

3) Steps 1-2 are repeated until all M units have been added. 

4) The Rosenbluth weight 

M 

= (4) 

i=l 

is calculated. This attempted move is accepted with a probability 

acc{o -^n)= min[l, W^^/iy . (5) 

The quantity W^""* is the Rosenbluth weight for the reverse move and is calculated as in 
steps 2-4, but with k — 1 random orientations and one orientation that is equal to the 
original geometry for each rigid unit. 

A type Ila move is very similar to a type I move. In this case, the direction of regrowth 
is chosen randomly. Then the first backbone unit to be regrown is chosen. The M rigid 
units to be regrown are removed and added back sequentially, as in the type I move. The 
rigid units in this case are either A-units, B-units with the side chain rigidly attached, C- 
units, or D- units (see appendix A). An alternative procedure would be to regrow the side 
chain units as well, but this proved not to be efficient, due to frequent steric repulsions. 
The move is accepted with the probability given by Eq. (^. 

A type lib move is initiated by identifying the 4 rigid units on the backbone to be 
rotated. This is done randomly. The four rigid units are labeled in an amine to carboxy 
terminal fashion. The attached side groups are rigidly rotated with the backbone units. 

The rotation is carried out as follows (see Fig. 2): 

1) The driver angle 0o is changed by an amount 6(j)o, where — A0 < 6(f)o < A(f). This is 
done k' times with probabilities according to the internal potential, 

p™*(0o,) exp[-f3ul^\4>o,)] . (6) 

The function u}^^{(j)oj) is the internal energy associated with this torsional angle. Only those 
values of 0o that lead to valid solutions for the modified torsional angles are considered. In 
the general case there will be a distinct 0i for each solution arising from the new value of 
00- Define fc^") to be the number of 0o-0i pairs. If = 0, the move is rejected. 

2) A 00-01 pair is picked with probability 

pTi<Po„ 01,) = exp[-/3M-*(0o,, 0i,)]/W^^") , (7) 
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where 

l^W = ^exp[-/?«-*(0o,,0ii)] . (8) 
i=i 

The function UQ^^{(j)oj,(f)ij) is the part of the external energy that couples this part of the 
backbone to the rest of the molecule. The value J*-"^ of the Jacobian is calculated for the 
new, chosen configuration (as detailed in Appendix B). 

3) The reverse move is considered. That is, a rotation about the new, chosen 0o-0i pair 
is considered, k' — 1 random values 50o are chosen. The original value of 0o is assigned to 
the fc'th value. This move results in k^°^ solutions for 0i. fc^"^ is always greater than zero, 
since the original configuration exists. (Special care is taken to ensure that the original 
configuration is found by the root finding procedure.) The Rosenbluth weight is assigned 
to W^°\ The value J^°^ of the Jacobian is also calculated for the original configuration. 

This attempted move is accepted with a probability 



accio n) = mm 



[1^ jW|y(n)/j(o)p^(o)| _ ^g^ 



Splitting the energy into internal and external parts is rather arbitrary. There are 
some constraints imposed, however, by the requirement that the normalization constants 



for Eqs. (|I]) and (|g) be independent of chain conformation |^. We assume for simplicity 
that = 0. One other natural choice, however, would set the internal part equal to the 
torsional terms in Hintra and set the external part equal to the rest of H . 

For any Monte Carlo scheme to properly sample the Boltzmann probability distribution, 
detailed balance must be satisfied. Refs. and prove that detailed balance is satisfied 
for the above scheme. 



3 Application to Polyglycine 

In this section we present the result of applying this configurational bias Monte Carlo 
method to two simple peptides, polyglycine Ge and constrained polyglycine CGeC. 

Figure 3 shows the energy of linear polyglycine as a function of Monte Carlo steps. This 
run took roughly 3 hours on a Silicon Graphics Indigo^. In Fig. 4 we show the end-to-end 
probability distribution for this system. Gaining this degree of convergence took a one-day 
run. 

Figure 5 illustrates the energy of the cyclic polyglycine as a function of Monte Carlo 
steps. This run took roughly 6 hours. Figure 6 provides a histogram of the number 
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of solutions found for each attempted concerted rotation. In rare cases the root finding 
procedure failed to find all the roots. In the construction of this plot, we rounded n^"^ 
up when it was odd. Figure 7 shows the histogram for the C/3SSC/3 dihedral angle, with 
the statistics taken from a run six times as long as that illustrated in Fig. 5. To give 
a feel for the barrier to rotation about this angle, we show in Fig. 8 the potential of 
mean force. This potential was determined by umbrella sampling This curve took 

two orders of magnitude longer to determine than did the probability distribution in Fig. 
7. The potential of mean force is contrasted with the energy associated purely with the 
C^SSC/3 torsional terms. Finally, Fig. 9 shows the result of classifying the configurations 
produced by the method into distinct stable conformations. Fuzzy clustering ||3^ was used 
to determine the dominant conformations, with the result that there are only two or three 
distinct conformations within this limited simulation run. The simulation run depicted in 
Figs. 5 and 9 took approximately 8 hours on a Silicon Graphics Indigo^. 

4 Discussion 

We see that with a very modest computational effort, we can achieve equilibrated results 
for linear peptides. With somewhat more effort, we can achieve equilibration for cyclic 
peptides. 

As expected, we find that the linear peptide Gq is relatively unstructured in solution. 
There is a common crumpled state, but there is also a significant population of the extended 
state. The constraint of the disulfide bond in CGeC, in contrast, forces that molecule to 
adopt a limited number of molecular conformations. For the fairly short runs illustrated 
in Figs. 5,6,7 and 9, we find only three dominant conformations. The first conformation is 
associated with the C/3SSC/3 torsional angle of 290°, whereas the other two are associated 
with angles of 88° and 98°. The first of these conformations is very tight, with 0.7 A 
fluctuations about the mean for all atoms in the molecule. The other two are somewhat 
looser, with roughly 1.2 A fluctuations. We see from Fig. 9 that even in this short run the 
method revisits previous conformations. In the limit of a long simulation, the time spent in 
each conformation would, of course, be proportional to the exponential of the free energy 
of the conformation. 

If CGeC were achiral, the potential of mean force in Fig. 8 would be symmetric about 
0° and 180°. Since the C^ carbons in the cystine residues are, in fact, chiral, the potential 
of mean force is not required to be symmetric. The asymmetry seen in Fig. 8 results from 
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the mean, chiral force of the rest of the molecule on the C/3SSC/3 torsion. In fact, the 
AMBER forcefield takes this chirality into account by reducing the symmetry of the 
carbon in cysteine. We have used this geometry The barrier at 0° is due to a high 
steric repulsion between the hydrogens on the carbons adjacent to the disulfide bond. 
This barrier is substantially higher than the barrier at 180°. 

From Fig. 8, we see that there is a very significant free energy barrier to rotation about 
the C/jSSC/j torsional angle. This figure was not constructed from a standard simula- 
tion run, but by the specialized procedure of umbrella sampling. It is clear from Fig. 7, 
however, that the present method is able to overcome this barrier and to properly sam- 
ple the relevant conformations even in a relatively short simulation. Any method such as 
molecular dynamics or standard Monte Carlo that makes only small, local changes to the 
configuration would never cross this barrier in a simulation of reasonable length. High 
temperature dynamics can allow systems to cross high barriers, but can not perform the 
requisite Boltzmann sampling to predict the physiologically relevant conformations. Only 
a biased method that makes fairly large geometrical changes is capable of dealing with 
such barriers in an automatic way, without resort to special techniques such as umbrella 
sampling. Furthermore, the ability to perform umbrella sampling has as a prerequisite 
the detailed knowledge of the important conformations and the paths between them. In 
our specific case, we find our method to be two orders of magnitude more efficient than 
umbrella sampling. 

5 Conclusion 

We have presented a Monte Carlo method capable of sampling the relevant room- or body- 
temperature configurations of linear and cyclic peptides. This method allows the study of 
peptides important in biological and technological settings. Our sampling of the disulfide 
dihedral angle in a prototypical cyclic peptide indicates that the method can explore widely 
separated regions of conformation space according to the proper Boltzmann distribution, 
even if the barriers between the regions are quite large. Previous simulation methods either 
fail to sample the proper thermal distribution or are vastly more computationally intensive 
and require detailed knowledge of the thermally accessible regions. The method can be 
extended to allow incorporation of explicit water molecules. The method can be extended 
to force fields with flexible bonds and angles. These extensions are subjects for future work. 
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Appendix A: Rigid Units 

As described, the algorithm assumes that bond lengths and angles are fixed. The only 
degrees of freedom, therefore, are torsional angles. Due to the extremely high force constant 
for rotation about a tt bond, even some torsional angles are fixed as well. An entire collection 
of atoms that is rigid is called a rigid unit. Such a unit has an incoming bond as well as 
several possible outgoing bonds. There are four backbone rigid units. Unit A is the starting 
NH3 group. Unit D is the terminal C00~ group. Unit B is the CqH group. Unit C is the 
CONH amide bond group. 

The residues are connected to the backbone by outgoing bonds from the B units. Table 
1 lists the decomposition of the amino acid side chains into rigid units. Typical rigid 
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Table 1: The rigid units in peptide side groups. 



Side CrrouD 


Risrid Units 


Glycine 


H 






A rrriTi iTip 


CHo CHo CHo CNqHp: + 






A G;"nP TP (TITIP 
JTi-iDlJCXiL CXpyllly^ 


CHo CONHo 


r^VGlf I P IITIP 


CHo SfHl 


frllltf^TTlJ^tP 


( ( H o ( 1 H o ( j( lo 

^■'-■'-z? v7vy2 


It n tI" mi Tl P 

VJ 1 LI U CUlllllC 


CHo CHo CONHo 


Histidine 






1 1 1 1 1 1 1 1 r» 1 1 1 1 r» 1 1 1 1 r> 

\^ii2) ^^^3) ^^^3 


T iPllPlTlP 




T (V<^ i Ti p 


CHo CHo CHo CHo NHq+ 

^^^2-) ^^^2-) ^^^2-) ^^^2-) ^^^^3 


Methionine 




P ri PTi \n Pi m m n p 




r roune 


JjaCKDOne vjroups. <^Qrl<^rl2>^rl2v>n.2) J^'l) v>W 


Serine 


CH2, OH 


Threonine 


CH, CH3, OH 


Tryptophan 


CH2, CgNHg 


Vahne 


CH, CH3, CH3 


Tyrosine 


CH2, C6H4, OH 



units are the CH2, CN3, CO2, and aromatic ring groups, which have substantial tt bonding 
character. 

Proline is a special case, technically an imino acid. The special nature is due to the 

cyclic bonding of the residue to the backbone. The rigid units in this amino acid are 
the CH„, CO, and N groups. Only trans isomers are allowed for the proline amide bond. 
Proline is treated in an approximate way: the C^-C^ fragment is kept rigid, the C^-N bond 
is broken, and the Cq.-N torsional barrier is increased. This approximation ignores the 
small fluctuations in the conflguration of the proline side-chain loop. 
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Appendix B: Concerted Rotation 



Since the molecules under consideration can be cyclic, a Monte Carlo move that preserves 
this constraint is required. The "concerted rotation" scheme used for alkanes p2| can 
be extended to allow rotation of the torsional angles in cyclic peptides. This appendix 
describes this extension. The reader is referred to Ref. for a fuller discussion of the 



original, restricted method. The method presented here allows for a fairly general molecular 
geometry. In particular, the method naturally accommodates the constraint of a planar 
amide bond. 

To formulate the method, we consider rotating about seven torsional angles, which will 
move the root positions of four rigid units, rotate up to three additional ones, and leave the 
rest of the peptide fixed. We define the root position of a rigid unit to be the Cq, position 
for a B unit, the C position for a C unit, the C position for a CH2 unit, and the S position 
for the S unit in cystine. If unit 5 is a C unit, however, is defined to be the N position 
of that unit. For each unit we define 6i to be the angle between the incoming and outgoing 
bonds. Thus, 6i = for a C unit, and 6i ^ 70.5° for all others. Figure 1 illustrates the 
geometry under consideration. 

The method leaves the positions of units i < or z > 5 fixed. The torsion (/)o is 
changed by an amount 6(j)Q. The values of 0j, 1 < i < 6, are then determined so that only 
the positions of units 1 < i < 4 are changed. 

The method required several definitions to present the solution for the new torsional 
angles. Vectors are defined which are the difference in position between unit i and unit 
i — 1, as seen in the coordinate system of unit i: 



Si) _ 

■i i—1 



(10) 



The coordinate system of i is such that the incoming bond is along the x direction. Thus 
Ij = /jX if atom and rj_i are directly bonded and has x- and y-components otherwise. 
We now define a rotation matrix that transforms from the coordinate system of unit i + 1 
to unit i 



I 



cos Qi 



sin Qi cos ( 



sin 



cos Qi cos ( 







sm I 



sin Qi sin 0j — cos di sin 0j — cos 0i y 
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The positions of the units in the frame of unit 1 are, thus, given by 



ll 

ll + T1I2 

Il + Ti(l2 + T2l3) 



r^^^ = Il + Ti(l2 + T2(l3 + T3l4)) . 



(12) 



We further define the matrix that converts from the frame of reference of unit 1 to the 
laboratory reference frame 



T'i«* = [cos^/^I + nn' (1 - cosip) + Msin?/;]A , 



(13) 



where 



M 



^ —Uz riy ^ 



Uz —n^ 



\ -riy rix 



(14) 



and 



n 



COS'0 = 

sin'^ = 



X X r 
|x X r| 
r • X 



|r| 

|r X x| 



(15) 



where x is a laboratory unit vector along the x direction, and r is the axis of the bond 
coming into unit 1. The matrix A is a rotation about x and is defined so that Ali = Ar: 

/ 1 ^ 



A = 



c -s 
s c 



(16) 



where 



= {hyAry + luAr,)/{Arl + Arl) 
= i-hAry + hyAr,)/{Arl + Ar^) 



(17) 
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Here Ar = A[T5^*^]^-'^(ri — ro) if unit is a C unit; otherwise Ar = li. 

The method proceeds by solving for 0j, 2 < i < 6, analytically in terms of 0i. Then a 
nonlinear equation is solved numerically to determine which values of (pi, if any, are possible 
for the chosen value of 0o- 

We will work in the coordinate system of unit 1, after it has been rotated by the chosen 
00- We define 

t = r(^)-li = [Tf]~i(r5-ro)-li. (18) 
If 6*3 7^ and 6*5 7^ 0, the square distance between unit 3 and unit 5 is known and equal to 



Qi = {hx cos 9/1 — liy sin 6^4 + /s^,)^ + {1^^ sin 6^ + l^y cos 6^ + / 



(19) 



But this distance can also be written as 



TrH - 1 



2 • 



(20) 



Equating these two results, two values of 02 are possible 



arcsinci — aicianXy/ Xz — H{xz) 
TT — arcsinci — arctanxy/x^ — H{ 



(21) 



with 



H(x) 




The constant Ci is given by 



Cl 



'' qj-X^-ll+2Xa:{cOS fe^Sa: +sin 02^31; ) n -J. f] Q _/ fi 

-2(sm02i3.-cose2«3y)(a;2+x2)l/2 ) t/3 ^ U, t/5 ^ U 

hx+Ux+kx COS d4-Xa: COS 02 Q _ Q ^ ^ Q 

(r5-r2)-(r6-r5)//6-<5z-U3: cos 6>4-3;3:(cos 62hx+sin62hy) 
(sine2'3a:-cos 92l'iy){xl+xl)^/'^ ' 

<3 ^ cos 6'4 -Xj; (cos 02 ^3 a: +sin 02 ^3 1; ) Q _ c\ o _ n 
(sine2hx-cose2hy){xl+x1)^/-^ ' - U, ^5 - U 



^3 7^0,^5 = 



(22) 



(23) 



where x is given by Eq. (|20|) if 6^ ^ 0, and x = T^\T^f]-\vQ - r^)/^ if 65 = 0. Clearly 
for there to be a solution |ci| < 1. The last three equations for ci were determined by 
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conditions similar to equating Eqs. (0) and (^). For 6*3 = 0,9^ 0, the x-component of 



,(3) 



Fg is known to be equal to ^.j. + I5 cos 6*4. For 63 ^ 0,65 = 0, the x-component of 



■5 

rl — ry is known to be equal to l^^ + ^42:Cos6'4. For 6^ = 0, 6^5 = 0, the angle between 
r3 — r2 and rg — is known to be equal to 64. 

To determine 03, two expressions for |r5 — r4p are again equated to determine 



and 



C2 



^^3 



II - 



ll + 2y^. (cos 6*3/4^ + sin6'3/4j 



-2(sin^3^4.-COSe3/4,)(l/2+y2)l/2 

arcsinc2 — arctan yy/yz — H{yz) 

71 — arcsinc2 — arctany^/y^^ — H{yz) , 



(24) 



(25) 



where y = ^{Ti — I2) — I3. Again, |c2| < 1 for there to be a solution. 
If 6*5 7^ 0, the value of 04 can be determined from 



,(1) 



rl') + T1T2T3T4I5 . 

q3 = T3-iT-^Tr^m-^ 



Defining 

the equations that define 04 are given by 



r5 - r4j 



(26) 
(27) 



(l3z 



cos 04 (sin 6^4/5^ — COs6'4/5j^) 

sin04(sin6'4/52, — cos6'4/5„) 



(28) 



This is a successful rotation if the position of rg is successfully predicted. That is, the 
equation 



T1T2T3T4T5I6 = [T'f]-\^e - r5 



(29) 



must be satisfied. We consider the x-component which implies 



,(1) M)\T 



T1T2T3T4X - (/e. cos 05 + hy sin ^5) = 0, ^5 7^ 



(r4 - r3) ■ (rg - rg) - ^4/6 cos 6*4 = 0, 6*3 7^ 0, = 



(30) 



Irg - r4| - 



1/2 



{k. + k.f + k; =0, ^3 = 0,^5 = 



must be satisfied if the rotation is successful. The equations for the case = clearly 
express the geometric conditions required for a successful rotation. 
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Eq. (pop is the nonlinear equation for 0i that must be solved. The equation depends 
only on 0i because 02, 03, and 04 are determined by Eqs. (pi)), (p5|), and (^) in terms of 
01. This equation has between zero and four values for each value of 0i, however, due to 
the multiple root character of Eqs. ( pi]) and (pS]). Equation (|30| ) is solved by searching the 
region — vr < < vr for zero crossings. The search is in increments of ~ 0.04°. These roots 
are then refined by a bisection method. There is always an even number of roots, due to 
the periodic nature of Eq. |30| . 

The root positions, rj, are enough to determine the position and orientation of the seven 
rigid units that are modified by the concerted rotation. Rigid unit is translated so that its 
root position is at tq. It is oriented so that its incoming bond vector is along the outgoing 
bond vector of rigid unit —1. It is then rotated so that its outgoing bond vector ends at 
ri. This process is repeated sequentially for rigid units 1 to 6. 

Repeated application of the concerted rotation leads to a slightly imperfect structure, 
due to numerical precision errors. In a practical application, the geometry would be restored 
to an ideal state by application of the SHAKE |^ or Random Tweak algorithm ^ . 



The transformation from 0j, < z < 6, to the new solution which is constrained to 
change only r^, 1 < i < 4, actually implies a change in volume element in torsional angle 
space. This change in volume element is the reason for the appearance of the Jacobian in 
the acceptance probability. The Jacobian of the transformation for alkanes is calculated 



in Ref. It is slightly different here since root position rs is not necessarily the head 



position. The Jacobian is given by 

J = l/|detB| , (31) 

where the 5x5 matrix Bij is given by the ith component of x (rs — hj) for z < 3 and 
by the {i — 3)th component of Uj x (rg — r5)/|r6 — rsl for z = 4, 5. Here hj = except that 
hs is the head position even if 6*5 = 0, and Uj is the incoming unit bond vector for unit i. 



Figure Captions 

Figure 1. The type I move applied to the serine side chain. 

Figure 2. The type lib move is illustrated for the case where unit is (a) a B-unit and 
(b) a C-unit. In each case, the original geometry and the four possible new geometries for 
the chosen driver angle are shown. In case (a), one of the new geometries is very different 
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from the original and the other three new ones. The move is shown for a hnear peptide, 
although it is used only on cyclic peptides. 

Figure 3. The energy of Gq as a function of Monte Carlo steps. Note the rapid 
equilibration. 

Figure 4. The probabihty distribution for the end-to-end distance for Gq. The distance is 
between the terminal Cq, groups. 

Figure 5. The energy of CGeC as a function of Monte Carlo steps. Note the rapid 
equilibration. 

Figure 6. The number of new solutions found for each attempted concerted rotation for 
CGeC. 

Figure 7. The observed probability distribution for the C^SSC^ torsional angle in CGqC 
is shown. 

Figure 8. The potential of mean force calculated by umbrella sampling for the C/3SSC/3 
torsional angle in CGgC (dashed line). The potential of mean force implied by Fig. 7 is 
indicated by the solid line. Also shown is the bare torsional energy contribution for this 
rotation (dotted line). 

Figure 9. Shown are the occupation numbers of the configuration in each of the three 
dominant conformations as a function of Monte Carlo steps (a). Also shown is the 
all-atom root-mean-squarc displacement of the configuration from each of the three 
dominant conformations (b). The curves for conformation 1 are solid, those for 2 are 
dashed, and those for 3 are short-dashed. 
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